function res = GetMainField(obj, varargin)
    import MeasurementComponent.Signal.*
    L = log4m.getLogger;
    
    p=inputParser;
    p.addParameter('resonant_frequency', 100e3, @(x) x>0);
    p.addParameter('driving_amplitude', 0.1, @(x) x>0);
    p.addParameter('guessField', 6.0, @(x) x>0);
    p.addParameter('fittingWidth', 1.0, @(x) x>0);
    p.addParameter('bz_min', 0.0, @(x) x>=0);
    p.addParameter('bz_max', 7.0, @(x) x>=0);
    p.addParameter('npoint', 100, @(x) x>0);
    p.addParameter('average', 100, @(x) x>0);
    p.addParameter('order', 4, @(x) ismember(x, 1:8));
    p.addParameter('bandwidth', 100, @(x) x>0);
    p.addParameter('isReport', false, @islogical);
    p.addParameter('description', 'normal', @ischar);
    p.parse(varargin{:});
    
    rbSignal = obj.setRFMagnetometer('field', p.Results.guessField, 'frequency', p.Results.resonant_frequency, 'amplitude', p.Results.driving_amplitude, ...
                                     'signal_oscillator', 1, 'signal_channel', 2, 'order', p.Results.order, 'bandwidth', p.Results.bandwidth);

    rbSignal.subscribe_signal(ziDemodTime.R_avg, [2 1 1]);
    rbSignal.subscribe_signal(ziDemodTime.Theta_avg, [2 1 2]);
    data = rbSignal.getParameterDependence(obj.ziSlave, ziSweepNode.AuxOut3_Offset, p.Results.bz_min, p.Results.bz_max, 'npoint', p.Results.npoint, 'averaging_sample', p.Results.average);
        
    [bzList, signalList] = data.getSweeperData(1);
    fitting = fitLorentzian(bzList, signalList, [p.Results.guessField-p.Results.fittingWidth, p.Results.guessField+p.Results.fittingWidth]);
    
    bz0 = fitting.bz0;
    obj.bz0.setDC(bz0); pause(0.5);
    L.info('ParametricMagnetometer::GetMainField', sprintf('Main field Control voltage is set to [%5.3f] V', bz0));
    phase0 = rbSignal.autoZeroPhase();
    L.info('ParametricMagnetometer::GetMainField', sprintf('Demod Phase is set to [%5.3f]', phase0));

    obj.by.disable();

    res.timestamp = datetime('now');
    res.RFSweepSignal = data;
    res.fitting = fitting;
    res.f0 = p.Results.resonant_frequency;
    res.bz0 = bz0;
    res.phase0 = phase0;
    res.width = fitting.width;
    res.hight = fitting.hight;
    obj.MainField = res;
    
    fig = obj.PlotMainField();
    %%
    if p.Results.isReport
        getReporter;
        rpt.NewRecord(['RFResonance_' p.Results.description]);
        rpt.AppendFigure(fig, sprintf('RFResonance %s', p.Results.description));
    end
    %%
    getSession;
    sess.addData('RFResonance', res, p.Results.description);
end

function fitting = fitLorentzian(bList, signalList, range)
    x = bList; 
    [maxY, pkIdx] = max(signalList);
    y2 = (signalList./maxY).^2;
    [xData, yData] = prepareCurveData( x, y2 );
    
    center = bList(pkIdx);
    range_width = range(2) - range(1);
    excludeIdx = find( x<range(1) | x>range(2));

    ft = fittype( 'a*g*g/((x-x0)^2+g*g)', 'independent', 'x', 'dependent', 'y' );
    excludedPoints = excludedata( xData, yData, 'Indices', excludeIdx);
    opts = fitoptions( 'Method', 'NonlinearLeastSquares' );
    opts.Display = 'Off';
    opts.Lower =      [0.5               0 range(1)];
    opts.StartPoint = [1.0 0.5*range_width center];
    opts.Upper =      [2.0 1.0*range_width range(2)];
    opts.Exclude = excludedPoints;

    [fitting.result, fitting.gof] = fit( xData, yData, ft, opts );
    fitting.bz0 = fitting.result.x0;
    fitting.width = fitting.result.g;
    fitting.hight = sqrt(fitting.result.a)*maxY;
    fitting.range = range;
    fitting.norm = maxY;
end
